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CONDITIONS FOR THE OCCURRENCE OF MEAN-MOTION 
RESONANCES IN A LOW MASS PLANETARY SYSTEM 

J.C.B. Papaloizou 1 and E. Szuszkiewicz 2 



Abstract. The dynamical interactions that occur in newly formed plan- 
etary systems may reflect the conditions occurring in the protoplane- 
tary disk out of which they formed. With this in mind, we explore the 
attainment and maintenance of orbital resonances by migrating planets 
in the terrestrial mass range. Migration time scales varying between 
~ 10 6 yr and ~ 10 3 yr are considered. In the former case, for which 
the migration time is comparable to the lifetime of the protoplane- 
tary gas disk, a 2:1 resonance may be formed. In the latter, relatively 
rapid migration regime commensurabilities of high degree such as 8:7 
or 11:10 may be formed. However, in any one large-scale migration 
several different commensurabilities may be formed sequentially, each 
being associated with significant orbital evolution. We also use a simple 
analytic theory to develop conditions for first order commensurabilities 
to be formed. These depend on the degree of the commensurability, 
the imposed migration and circularization rates, and the planet mass 
ratios. These conditions are found to be consistent with the results of 
our simulations. 



1 Introduction 

The increasing number of extrasolar multi-planet systems have stimulated stud- 
ies of their origin, evolution and stability. An important feature in planetary 
system evolution is the occurrence of mean- motion resonances, which may relate 
to conditions at the time of or just after the process of formation. There are 
some well known examples of systems such as Gliese 876 \ Marcy et al. 2 001), HD 
82943 d Mayor et al. 2004[ ) and 55 Cancri (| McArthur et al. 2004P involving plan- 
ets with masses in the Jovian range. More recently a system of three superearths 
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orbiting HD40307 has been announced for which the period ratios are not strictly 
commensurable but are close to a 2:1 commensurability Q Mayor et al. 2 008). 

It is therefore important to establish the main features of the evolution of low 
mass planets embedded in a gaseous disc and in particular to determine the types 
of resonant configurations that might arise when a pair of such planets evolves 
together. The disc-planet interaction naturally produces orbital migration through 
the action of tidal torques (jGoldreich fc Tremaine 19801 |Lin fc Papaloizou 1986| ) 
which in turn may lead to an orbital resonance in a many planet system (eg. 
Sncllgrove, Papaloizou & Nelson 2001 Lee & Pcale 2002). This is because it is 



expected that two planets with different masses will migrate at different rates. 
This has the consequence that their period ratio will evolve with time. In the 
situation where the migration is such that the orbits converge, they tend to enter 
and become locked in a mean-motion resonance (eg. |Nelson fc Papaloizou 2002 
|Kley, Peitz fc Bryden 2004| and then subsequently migrate together. 
In the simplest case of nearly circular and coplanar orbits the resonances that are 
formed are the first-order resonances which occur at locations where the ratio of 
the two orbital periods can be expressed as (p + l)/p, with p being an integer. 

Papaloizou & Szuszkiewicz (2005) performed a recent analytic and numerical 
study of the formation of first order commensurabilities by a system of two planets 
in the earth mass range migrating in a laminar disc. Here we extend these studies 
using TV-body simulations to a larger range of migration rates and commensu- 
rabilities and compare the numerical work to analytically derived conditions for 
particular commensurabilities to form. We begin by deriving the analytic criteria 
and go on to present the results of simulations to which they are reconciled. We 
consider migration time scales ranging between ~ 10 3 yr and ~ 10 6 yr. For the 
longest time scales, 2:1 commensurabilities may be set up while for the shortest, 
commensurabilities as high as 11:10 have been found. Finally we summarize our 
results. 



2 Basic Equations 



The equations of motion for a system of N planets (i = 1, 2, 3, ....N) moving in 
a fixed plane, about a dominant central mass, under a general Hamiltonian H may 
be written in the form (sec e.g. Papaloizou 2003, Papaloizou & Szuszkiewicz 2005) 
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Here the orbital angular momentum of planet i which has reduced mass m, = 
rriioM / (M + mio), with mm being the actual mass, is Li and the orbital energy is 
Ei. For motion around a central point mass M we have 



Lt = 



m^GM.a^l - ef) and (2.5) 

GMim, 

Ei = (2.6) 

Here Mj = M + mio, a,i denotes the semi-major axis and ej the eccentricity of 
planet i. 

The mean longitude of planet i is A,; = rii(t — toi) + v&i, where rii = y GM^ /o| 
is its mean motion, with toi denoting its time of periastron passage and Wi the 
longitude of periastron. 

In this paper we consider N = 2 and arrange the two planets such that i = 1 
denotes the outer planet and i = 2 denotes the inner planet. We remark that the 
above formalism does not incorporate disk tides. However, forces resulting from 
these can be added in separately (see below). 

2.1 Coordinate system 

We adopt Jacobi coordinates (jSinclair 19751 [Papaloiz ou fe Szuszkiewicz 2005| ) for 
which the radius vector r2, is measured from M and that of the outer planet, ri, 
is referred to the centre of mass of M and the inner planet i = 2. The required 
Hamiltonian correct to second order in the planetary masses is given by 

I, ,. 2 , ,. l2 , GMxmi GM 2 m 2 



H = +m 2 \r 2 \ ) 

Fil l r 2| 
,3 ""■ ( 2 -7) 



2 ki| \r 2 \ 

Gm\m 2 Gm\m 2 T\ ■ r 2 



Fia\ l r il 
Here Mi = M + mi , M 2 = M + m 2 and r 12 = r 2 — r! . 

The Hamiltonian may quite generally be expanded in a Fourier series involving 
linear combinations of the three angular differences Aj — w^, i — 1,2 and wi — vj 2 . 

Near a first order p + 1 : p resonance, we expect that both tfii = (p + l)Ai — 
p\ 2 — w 2 , and <j) 2 = (p + l)Ai — pA2 — tni, will be slowly varying. Following 
standard practice (see e.g. [Papaloiz ou 2003 1 and |Papaloizou fc Szuszkiewicz 2005) 
only terms in the Fourier expansion involving linear combinations of <pi and 4> 2 as 
argument are retained because only these are expected to lead to large long-term 
perturbations. In general there are an infinite number of such terms that need to 
be considered. However, in the limit of small eccentricities only terms that are first 
order in the eccentricities need to be retained. This approximation is valid when 
the circularization times are small enough to ensure that the eccentricities remain 
small. This situation is realized for realistic examples of low mass protoplanets 



migrating in protoplanetary discs (Papaloizou fc Szuszkiewicz 2005) 



Following the procedure described above and expanding to first order in the ec- 
centricities, the Hamiltonian may be written in the form H = Ei + E 2 + Hi 2 , 



4 



Title : will be set by the publisher 



where 



with 



H 



d = 



at 

d{b%{x)) 
dx 



(eiCi cos(<£ 2 ) - e 2 C 2 cos(</>i)) , 



+ (2p + l)6^ 2 (a;) - 4a;<5? I and 



(2.8) 



(2.9) 



C 2 = 



dx 



+ 2{p+l)b^'(x) 



(2.10) 



Here 6^ 2 (ar) denotes the usual Laplace coefficient (e.g. Brouwer & Clemence 1961) 
with the argument x — a 2 /ai and 5f denotes the Kronecker delta. We have also 
replaced M t by M. 



2.2 Behaviour near a resonance with disk tides incorporated 

The governing equations for motion near to the p + 1 : p resonance are to lowest 
order in the eccentricities 



I de\ 



Gm 2 C 1 



= ei sin 02 



tcl 



-.e 2 sin (pt 



2 dt aWGMci! 
1 de\ GmiC 2 
2~dt~~ ai y/GMa 2 

"&ni{p + l)Gm 2 
ax^/GMal 

3n 2 pGm 1 . 

= (Giei sin <p 2 - C 2 e 2 sin </>i) + 



(Ciei sin 2 - C 2 e 2 sin + 



n 2 = 



aiy/GMa.2 



t 
n 2 



t c i 

2*4 



3n 2 e 21 



tmig2 



t. 



c2 



(2.11) 
(2.12) 
(2.13) 
(2.14) 



Here the terms on the right-hand sides enclosed in square brackets give the contri- 
butions arising from disk tides. These have been discussed elsewhere (see e.g. 
|Papaloizou fc S zuszkicwicz 2005}. The circularization and migration times for 
planet i are t ci and t migi respectively. Note that the migration time is defined 
here as the time for the mean mo tion t o cha nge by a factor of e as a result of disk 
torques. The terms tx ef/t c i in (2.13) and (2.14) account for the orbital energy 



dissipation occurring as a result of circularization at the lowest order in a which 
it appears. 



The other terms on the right-hand sides of (2.11 1- (2. 14 1 are derived from the 



Hamiltonian system (2.1)-(2.4) using the Hamiltonian (2 



2.3 Migration maintaining resonance 

When the two planets migrate together maintaining resonance, we expect and find 
solutions for which e\ and e 2 are either actually nearly constant or more generally 
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very nearly constant in a time average sense while the ratio n\jni is maintained 
very close to p/(p + 1). When considering time averages, we here consider the 
average to be taken over a time long compared to the characteristic orbital period 
but short enough that the semi-major axes and tidal time scales may be consid- 
ered constant. Using angle brackets to denote such a time average and setting 
(eiei) = (e-ie-i) = in (2.11) and (2.12) and taking the result of dividing the 



time averaged (2.13) by the time averaged (2.141 gives three equations from which 
(ei sin(0i)), £ = 1,2 may be determined together with a constraint on the mean 
square eccentricities and migration times that does not depend on fa. Proceeding 
in this way we obtain 



(eism d> 2 ) 



(e 2 smfa) 



ela-Ly/GMal 
GmiC 2 t C 2 



(2.15) 
(2.16) 



together with the constraint 



e\ e 2 rri2n\a\ 



tc 



t C 2 TOin 2 a2 



p 2 



tr 



migl 



t 



migl 



(2.17) 



where / = m2a\/((p + l)(m2ai + toici^)) and for ease of notation we have sup- 
pressed the angle brackets enclosing e 2 which is now implicitly assumed to be 
a time averaged value. We remark that (2.171 has already obtained from gen- 



eral considerations ( |Papaloizou fc Szuszkiewi cz 2005] ) . However, (2.151 and (2.161 
have not. 



2.4 A condition for a p+l:p commensurability to be maintained 



We m ay use (2.15 1 and (2.16 1 to express e~ in terms of (e^ sin fa) in the constraint 
(2.171. Using Schwartz's inequality and applying the condition (sin 2 fa) < 1 results 
in an inequality that must be satisfied if the resonance is to be maintained. This 
inequality can be written in the form 



p 2 n^m2 
(p+ 1) 2 M 



(1 - f)m 2 C(t cl miafCff c2 



7V 



Mai 



> 



1 



1 



migl 



migl 



(2.18) 



For the simple example where mi 3> m 2 is in a prescribed slowly shrinking circular 
orbit and controls the migration {t mig 2 S> t mig i) the relation (2.181 simplifies to 
the form 



M 2 



> 



"2 



3pa 1 n 1 n2t mlgl t C 2C 2 



(2.19) 



We further remark that when it is satisfied, the maximum (sin 



has to exceed 



the ratio of the right-hand side to the left-hand side of the inequality (2.18 1. 
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2.5 Resonance overlap 

Because it is found that both C\ and C2 increase with p, while / decreases with 



p, the inequality (2.18) indicates that for given planet masses the maintenance 
of resonances with larger values of p is favoured. However, the maintenance of 
resonances with large p may be prevented by resonance overlap and the onset of 
chaos. Resonance overlap occurs when the difference of the semi-major axes of 
the two planets is below a limit that, in the case of two equal mass planets, has a 



half- width given by Gladman (1993) ] as 

with a and m p i ane t being the mass and semi-major axis of either planet respec- 
tively. Thus for a system consisting of two equal four-Earth-mass planets orbiting 
a central solar mass we expect resonance overlap for p > 8. Conversely we might 
expect isolated resonances in which systems of planets can be locked and migrate 
together if p < 8. But note that the existence of eccentricity damping may allow 



for somewhat larger values of p in some cases. In this context the inequality (2.18 I 
also suggests that resonances may be more easily maintained for lower circulariza- 
tion rates. However, this may be nullified for large p by the tendency for larger 
eccentricities to lead to greater instability. Note also that higher order commcn- 
surabilities may also be generated in such cases and these are not covered by the 
theory described above. 



3 Numerical simulations 



In order to check the applicability of the inequality (2.18) we have performed nu- 
merical simulations of a pair of migrating planets using an iV-body code. The 
approach is the same as that used by many authors (e.g. Snellgrove, Papaloizou & 
Nelson 2001; |Nelson fc Papaloizou 20021 [L ee & Peale 20021 Kley, Peitz & Bryden 
2004). The reader is referred to these papers for details. In particular, orbital 
migration and eccentricity damping that are presumed to result from interac- 
tion with the protoplanetary disk are modelled through the addition of appropri- 
ate acceleration terms to the equations of motion. We considered planets with 
masses mi = 4M®, and rri2 — M® for a wide range of imposed migration and 
circularization rates. The planets were started in circular orbits with mi be- 
ing at 0.34aw and ni2 at 0.2au (note that these results can be scaled to apply 
to other initial radii by appropriate renormalization of space and time scales). 
Thus they begin just outside of the 2:1 resonance before commencing conver- 
gent migration. Low mass planets of the type we consider undergo type I mi- 
gration (see Papaloizou &; Terqucm 2006 and references therein). On account of 



dependence on the equation of state and non linear effects, the precise rate that 
should be employed is uncertain even in a laminar disk as adopted here (e.g. 
|Paardekooper fcPapaloizou 2008J . To survey a range of possibilities, the migra- 
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Table 1. Details of the simulations discussed in this paper. Column 1 gives the value 
of /tid- Column 2 gives the run time in years. Column 3 lists the commensurabilities 
for which sustained trapping and joint migration was noted. The occurrence of such 
episodes was characterized by eccentricity boosts. Column 4 firstly indicates the first 
order commensurability of lowest degree estimated to be possible by use of the inequality 
(2.181 and secondly the same estimate assuming the maximum allowed (sin 2 <f>i) < 0.1. 
All the simulation results are consistent with the inequality (2.181. 



/tid 


Run time yr. 


Commensurabilities 


Lowest degree expected 


1 


1278400 


2:1 3:2 


2:1, 2:1 


2 


638400 


2:1 3:2 


2:1, 2:1 


4 


638400 


2:1 3:2 


2:1, 2:1 


16 


492800 


2:1, 3:2, 4:3 


2:1, 2:1 


32 


345600 


2:1, 3:2, 4:3 


2:1, 2:1 


64 


123200 


3:2, 4:3, 5:4 


2:1,2:1 


128 


57000 


3:2, 4:3, 5:4, 6:5 


2:1, 3:2 


256 


31150 


4:3, 5:4, 6:5 


2:1, 3:2 


512 


15550 


4:3, 5:4, 6:5 


3:2, 5:4 


1024 


7875 


5:4, 6:5, 7:6 


4:3,7:6 


2048 


3906 


6:5, 7:6, 8:7 


5:4, 10:9 


4096 


1959 


11:10 


8:7, p > 12 




400000 800000 1 .2e+0S 1 .6e+06 
time yr. 




400000 800000 1 .2e+06 1 .6e+0E 
time yr. 



Fig. 1. The evolution of the semi-major axes (left panel) and eccentricities (right panel) 
of the two planets are shown as functions of time for / t id = 1. The planets migrate inwards 
while locked in commensurabilities. The early episode of high eccentricity corresponds to 
a 2:1 commensurability while the final one corresponds to a 3:2 commensurability. The 
inner planet has the larger eccentricity. 
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120000 160000 




40000 80000 120000 160000 

time yr. 



Fig. 2. As in Fig. [I] but for / t id = 64. In this case the planets fail to lock into the 2:1 
commcnsurability, showing only an eccentricity spike when the system passes through 
this commensurability. Subsequently the configuration becomes sequentially locked in 
the 3:2, 4:3 and 5:4 commensurabilities. 



0.4 | 1 1 1 1 1 0.035 




1000 2000 3000 4000 5000 1000 2000 3000 4000 5000 

time yr. time yr. 



Fig. 3. As in Fig. [T]but for /tid = 1024. In this case the planets fail to lock into the 2:1, 
3:2 and 4:3 commensurabilities showing only eccentricity spikes when the system passes 
through these commensurabilities. Subsequently the configuration becomes sequentially 
locked in the 5:4, 6:5 and 7:6 commensurabilities. 



tion and circularization rates we adopted were given by 



Papaloizou & Szuszkiewicz: MMRs in a low mass planetary system 



9 



0.3 | 1 1 1 1 1 0.025 




200 400 600 800 1000 200 400 600 800 1000 

time yr. time yr. 

Fig. 4. As in Fig. [T] but for /yd = 4096. In this case the planets fail to lock into any 
commensurability until p — 10, corresponding to the 11:10 commensurability. Many 
eccentricity spikes occur as the system passes through lower degree commensurabilities. 



and 

_ 2.5 x 10 4 M e 
/tid mi 



yr, (3.2) 



respectively where / t id is a scaling factor. We have adopted 1 < / tic j < 4096, 
giving migration time scales that range between ~ 10 6 yr and ~ 10 3 yr. In this 
way timescales that are both comparable to the gas disk lifetime, and significantly 
shorter are considered. The ratio t C i/t m i g i has been chosen to be comparable to the 
disk aspect ratio squared for an aspect ratio ~ 0.05 as expected from theoretical 
work (see e.g. |Papaloizou fc Larwood 2000| 

We summarize the parameters and outcomes of the simulations in Table [l] In par- 
ticular the commensurabilities for which sustained trapping and joint migration 
were observed are noted. As expected, it is found that lower degree commen- 
surabilities are formed for smaller / t id and then become unsustainable if / t id is 
increased beyond a critical level. In general each simulation is found to involve the 
setting up of a sequence of commensurabilities each of which is disrupted after a 
significant contraction of the planetary orbits. 

To illustrate this, the evolution of the semi-major axes and eccentricities of the 
two planets for /tid — 1 are plotted in Fig. [I] In this case a 2:1 commensurability 
is first formed which is then disrupted after the planets have migrated inwards 
together contracting their orbits by ~ 30%. This resonance is disrupted and a 3:2 
commensurability is then formed. The evolution for /tid = 64 is illustrated in Fig. 
[2j In this case the planets pass through the 2:1 commensurability. Subsequently 
the configuration sequentially forms 3:2, 4:3 and 5:4 commensurabilities. Each 
of these are associated with significant contractions of the orbits. The case with 



10 



Title : will be set by the publisher 



/tid = 1024 is illustrated in Fig. [3] In this case the system passes through the 2:1, 
3:2 and 4:3 commcnsurabilities. Subsequently 5:4, 6:5 and 7:6 commensurabilitics 
are formed each of which is associated with a significant orbit contraction. The 
results for / ti( j = 4096 are shown in Fig. |4j No lasting commensurability is formed 
until p = 10. Again this commensurability survives while the orbits undergo a large 
contraction. This commensurability is of somewhat higher degree than would be 
expected from |Gladman (1993) ~| It is likely that it can survive on account of the 
eccentricity damping employed here. 

In Table [l] we also indicate the first order commensurability of lowest degree 
that is estimated to be possible by use of the inequality (2.181 together with an 
estimate made assuming the maximum allowed (sin 2 (/>,;) < 0.1. All the simulation 
results are consistent with the inequality (2.181. If it is used directly, the lowest 
degree commensurability that is possible to form is, as expected, somewhat un- 
derestimated. However, if we assume the maximum allowed (sin 2 4n) < 0.1, Table 
[l] indicates that the lowest degree commensurability that is possible to form is 
somewhat overestimated for the larger values of /tid- 



4 Conclusion 

We have studied the development of orbital resonances by migrating planets in 
the terrestrial mass range for migration time scales varying between ~ 10 6 yr 
and ~ 10 3 yr. In a typical simulation, a sequence of resonances occurred, each 
of which could survive significant orbital evolution before being lost. For the 
slowest migration rates considered, a 2:1 commensurability was able to form. At 
the fastest rates only commensurabilities with large p persisted. We also found 
analytic conditions for first order commensurabilities to be formed. These were 
consistent with our numerical simulations. 
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